• cholera is waterborne disease

• the Broad Street pump is the source of the outbreak

snowMap(); addKernelDensity(obs.unit = "fatality")

• does the map support his claims?

• what would waterborne transmission look like?
• what would airborne transmission (miasma) look like?
• does the map help us to distinguish the two?

"The inner dotted line on the map shews [sic] the various points which have been found by careful measurement to be at an equal distance by the nearest road from the pump in Broad Street and the surrounding pumps …"

"… it will be observed that the deaths either very much diminish, or cease altogether, at every point where it becomes decidedly nearer to send to another pump than to the one in Broad Street".

• Report On The Cholera Outbreak In The Parish Of St. James, Westminster, During The Autumn Of 1854, p. 109.

• if cholera is waterborne & …
• if pumps are the source of water,
• then the extent of outbreak is limited by a pump's reach

• the set of homes defined by proximity to a pump
• where TO find cases (inside)
• where NOT to find cases (outside)

• proximity matters
• fetching water from a well is a chore
• so choose the nearest pump
• (rival explanations)

• "choose the nearest" implies choice
• other pumps shape a pump's neighborhood

• hand drawn, eyeballed quality

• what about other neighborhoods?
• what about hypothetical/counterfactual neighborhoods?

pumps

##    id                    street         x         y
## 1   1              Market Place  8.651201 17.891600
## 2   2        Adam and Eve Court 10.984780 18.517851
## 3   3            Berners Street 13.378190 17.394541
## 4   4             Newman Street 14.879830 17.809919
## 5   5          Marlborough Mews  8.694768 14.905470
## 6   6 Little Marlborough Street  8.864416 12.753540
## 7   7              Broad Street 12.571360 11.727170
## 8   8            Warwick Street 10.660970  7.428647
## 9   9             Bridle Street 13.521460  7.958250
## 10 10             Rupert Street 16.434891  9.252130
## 11 11               Dean Street 18.914391  9.737819
## 12 12          Tichborne Street 16.005110  5.046838
## 13 13               Vigo Street  8.999440  5.101023

deldir::deldir(pumps[, c("x", "y")])

plot(deldir::deldir(pumps[, c("x", "y")], suppressMsge = TRUE))

sp::point.in.polygon(cholera::fatalities.address$x,
  cholera::fatalities.address$y, cell$x, cell$y)

plot(neighborhoodVoronoi())

plot(neighborhoodVoronoi())

plot(neighborhoodVoronoi(-6))

plot(neighborhoodVoronoi(-7))

plot(neighborhoodVoronoi(6:7))

• fast (uses only location of pumps)

• doesn't quite match Snow's annotation

plot(neighborhoodVoronoi()); addSnow()

• uses only location of pumps
• Euclidean distance (Snow implies walking distance)
• network of roads and streets don't matter
• walk through walls

plot(neighborhoodEuclidean(case.set = "observed")); addVoronoi()

• compute and plot the path between case and pump …

• do so along roads and streets of on the map …

• rinse and repeat -> pump neighborhood

• Dodson and Tobler (1992)
  • GIS exercise
• bars (cases) and pumps are points
  • x-y coordinate
• streets are approximated by straight line segments
  • x-y coordinates of segment endpoints

• relationship between case and road ("address")
• relationship between cases in stack

orthogonal.slope <- -1 / segment.slope

orthogonal.intercept <- case$y - orthogonal.slope * case$x

x.proj <- (orthogonal.intercept - segment.intercept) /
          (segment.slope - orthogonal.slope)

y.proj <- segment.slope * x.proj + segment.intercept

• label/group cases as a stack

stats::hclust()
stats::cutree()

• roads -> edges
• street intersections -> nodes (vertices)

plot(neighborhoodData(embed = FALSE))

plot(neighborhoodData(embed = TRUE))

neighborhoodData <- function(vestry = FALSE, case.set = "observed",
  embed = TRUE, embed.landmarks = FALSE) {

  node.data <- nodeData(embed = embed, embed.landmarks = embed.landmarks)
  nodes <- node.data$nodes
  edges <- node.data$edges
  g <- node.data$g
  nodes.pump <- nodes[nodes$pump != 0, ]
  nodes.pump <- nodes.pump[order(nodes.pump$pump), c("pump", "node")]
  out <- list(g = g, nodes = nodes, edges = edges, nodes.pump = nodes.pump)
  class(out) <- "neighborhood_data"
  out
}

plot.neighborhood_data <- function(x, ...) {
  plot(x$g, vertex.label = NA, vertex.size = 2, ...)
}

plot(neighborhoodData())

• find nodal path between case and pump
• many possible paths
• use the "shortest" path
• translate nodal path -> walking path

# breadth first algorithm
igraph::shortest_paths(g, ego.node, alter.node)

walkingPath(447, weighted = FALSE)$path

##  [1] "14.5236060046498-11.3601587685374"
##  [2] "14.6879997-11.0889997"            
##  [3] "14.8257446-10.850667"             
##  [4] "14.9709997-10.5930004"            
##  [5] "15.1669998-10.3299999"            
##  [6] "15.5965776-10.6166334"            
##  [7] "16.0209999-9.8280001"             
##  [8] "16.3169994-9.4689999"             
##  [9] "16.9459991-10.0389996"            
## [10] "17.0279999-9.9549999"             
## [11] "17.2169991-9.7259998"             
## [12] "16.552-9.0880003"                 
## [13] "16.4352026365887-9.25235128027873"

plot(walkingPath(447, weighted = FALSE), zoom = FALSE, unit.posts = NULL)

# Dijkstra algorithm
igraph::shortest_paths(g, ego.node, alter.node, weights = edges$d)

plot(walkingPath(447, weighted = TRUE), zoom = FALSE, unit.posts = NULL)

plot(neighborhoodWalking())

plot(neighborhoodWalking())

plot(neighborhoodWalking(-6))

plot(neighborhoodWalking(-7))

plot(neighborhoodWalking(6:7))

• segment endpoints & …
• endpoints of split segments
• cut points via base::rle()

plot(neighborhoodWalking(case.set = "expected"), type = "road")

• observed paths -> expected areas
• 20K simulated cases (points)

sp::spsample(sp::Polygon(map.frame[, c("x", "y")]), n = 20000,
  type = "regular")

lapply(list, function)

mclapply(list, function)

head(roads[roads$name != "Map Frame", ])

##    street n         x        y id               name
## 26     13 2  9.682715 18.20924 25 Castle Street East
## 25     13 2 10.795917 18.56957 26 Castle Street East
## 31     15 2  8.657000 18.08900 30       Market Place
## 30     15 2  8.533000 18.52400 31       Market Place
## 33     16 2  5.097822 18.05334 32    Margaret Street
## 32     16 2  6.478946 18.47905 33    Margaret Street

head(road.segments)

##    street   id               name       x1       y1        x2       y2
## 26     13 13-1 Castle Street East 9.682715 18.20924 10.795917 18.56957
## 31     15 15-1       Market Place 8.657000 18.08900  8.533000 18.52400
## 33     16 16-1    Margaret Street 5.097822 18.05334  6.478946 18.47905
## 34     17 17-1      Regent Street 4.992000 18.44100  5.097822 18.05334
## 37     18 18-1 Castle Street East 8.772448 17.91459  9.682715 18.20924
## 38     19 19-1       Market Place 8.772448 17.91459  8.657000 18.08900

road.segments <- parallel::mclapply(unique(rd$street), function(i) {
  dat <- rd[rd$street == i, ]
  names(dat)[names(dat) %in% c("x", "y")] <- c("x1", "y1")
  seg.data <- dat[-1, c("x1", "y1")]
  names(seg.data) <- c("x2", "y2")
  dat <- cbind(dat[-nrow(dat), ], seg.data)
  dat$id <- paste0(dat$street, "-", seq_len(nrow(dat)))
  dat
}, mc.cores = cores)

road.segments <- do.call(rbind, road.segments)

simulateFatalities(compute = TRUE, multi.core = FALSE)

simulateFatalities(compute = TRUE, multi.core = TRUE)

plot(neighborhoodWalking(case.set = "expected"), type = "area.points")

plot(neighborhoodWalking(case.set = "expected"), type = "area.polygons")

• vector graphics

streetNameLocator("marshall street", zoom = TRUE)
addNeighborhoodCases(type = "expected")

streetNameLocator("marshall street", zoom = TRUE)
addNeighborhoodWalking()